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A complete automated computer program was prepared and hand delivered to 
NASA, Langley. The program Included complete listings, a deck of cards, and 
several detailed examples, all In the form of actual computer outputs. 

Mr. Isaac Horowitz, in a visit to NASA, Langley, gave an outline of the 
program. Mr. Patrick Rosenbaum spent several days at NASA, Langley, advising 
and assisting NASA personnel in its use by working with them at the NASA 
computer. A separate outline of the program, including explanation of the 
various steps, references to the mathematical techniques, used, etc., is 
included herein. 

Present efforts are now in the improvement and refining of the more 
sophisticated parts of the program -- primarily, the automatic evaluation 
of a sub-optimal loop transmission, and the translation of time to frequency 
domain specifications- 
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AUTOMATIC DESIGN. PROGRAM FOR SINGLE -LOOP SYSTEMS WITH UNCERTAIN 
■ PLANT PARAMETERS (ADSUP) 


Introduction 


This program is divided into four parts: 

PART 1 . is the translation of time domain boundaries (T.D.B.) into frequency 
•domain bounds (F.D.B.). 

PART 2 computes, for a given plant,, the templates at various frequencies and 
stores them, 

PART 3 derives the optimum L (jw) (open-loop transfer function with excess 

• ■ 2 

of 2 poles over zero). which satisfies the bounds on L(joo) due to 

F.D.B. , 

\ 

S 

PART 4 derives the final L (jo)) (excess of 5 poles over zero). It then 

5 

derives the prefilter F(jo).) . Both , L and F are obtained as 
rational functions in the complex variable s . 



2 


I. PART 1 


Philosophy 

In the first step; all second order systems (two poles, no zeros) are 
found, which satisfy the T.D.B. Real poles and zeros are then added, which 
•continue to satisfy the T.D.B. Part 1 consists of the following components. 

A. STARTl 


This routine reads the T.D.B. data. 

NT : number of T.D.B. points. 

DELTA: increment (in seconds) in time . 

BMAX ; vector of NT .maximum tolerable values. 

BMIN : vector of NT minimum tolerable values. 

Note : If NT<24 linear interpolation is used in order to have at least 24 
points. A maximum of 100 points is allowed (NT^lOO) . 

B. TWO 


This routine finds all second order systems which satisfy the T.D.B. 

First,- KRANGE computes and w corresponding to 

Jiiiii n iHdx 


03 * < CO < (0 

n mm n n max 


such that Vt BMIN(t) 


03 


s(s^+ao^s-to^) 


-) < BMAX(t) 



3 


The vectors TMIN{’) and TMAX{») are calculated on a set of selected 
frequencies (see TABCS). TMIN and TMAX denote respectively the minimum 
and maximum permitted frequency domain values of the system transfer function. 
Then, for e [u , there is found the range (by calling KRANGE) 

of the damping denoted CMIN and c ^ denoted CMAX] , which satisfies 

in 1 n rnaX 

the T.D.B.; is incremented using step-size DW . At each step TMIN and 
TMAX are, if necessary, updated. 

Note ; WO is called the central frequency; If M=128 (number of points) (see 
TABCS) then W0^W(65) . In general ^'W0^W(^+ 1) . 


1) KRANGE (A.,B„N) 

) 

t 

Input : A,B are the first guess for the range, N is number of iterations 

allowed. 

Output: A,B is the range of variation. 

The method used here is systematic cut and try. 


2) MTEST (X) 


Let X represent or c : 

If KFLA6=1 (in the common block), then £=1 

T(s) = if e T.D.B. 

+ 2 - + 1 

i . X ■ 


and we test the transfer function 


iPRODUCIBiLrrY OF W 
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If KFLAG=2 then we -test if the step response of T(s) 
lies inside the time domain bounds. 


(OMEGA) 


2 OReM ^ 


3) KTEST (NQ,NP.Q,P,GAIN,WO) 


We test if GAIN*Q{S/W0)/P{S/W0) lies in T.D.B. Q is a polynom of degree 
NQ . P is a polynom of degree NP . The outputs are shown in Fig. 1.1. 



C.‘ STEP and dSTEP 


Those routines compute the step response of G*Q(S)/P{S) . 

' / 

1) -STEP (M,N,Q,P,G,DELTA) 

/ 

Q is a polynom of degree M . P is a. polynom of degree N 


Let the dynamical system be given by; 
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b = 


X = ApX + bu 


y = cx 


0 

G 


where: A„ = C 


-Pr 


••P. 




and let u be a step function. Then y(s) = i C(SI-A )~^b . H(S) = C(SI-A l^^b 

•J u u 


is the Laplace transform of the impulse response. 


5A. 


If . 6=DELTA , is the increment in time, let A=e ^ . The step response 
i s then : 

t A„(t-x) 

y(t) == / Ce P bdT 
0 

At. 

• y{t) = -CAp b + Ce P Ap'b 

A t 

= H(0) + Ce P Ap^b 


but 




b = 


-g/Pn . 


so finally: 


y(k6) = H(0) + (-g/p^)(CA'')^ 


( 1 .) 


where the subscript i refers to the first element of that vector. 

I 

The following notation is used: 

A <= EXP(DELTA*Ap) 

B <= -g/p^ 

PN <= P(N‘) 

F <= H(0;) = q/g/p^ 
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AA<= Ap 

STEP is the value at t=0 of the step response, 
2). DSTEP (N) 


N is the degree of P(*) . To calculate the step response at the time 
kx6 , we must call k times the function DSTEP , which computes (1). 

D. Description of the Package of Subroutines for Stable, Minimum Phase 
Rational Approximation to a Magnitude Squared Function (RFA) 

\ 

1) Philosophy 

The problem is to construct a rational function 
H(s) = g^ 

whose poles and zeros are in the LHP , such that 

|H(.jw)|'2 F (u) 

0 

where F (w) is a given magnitude squared function with a known excess of 
0 

poles over zeros. 

2) Theoretical' Preliminaries 

\ 

The heart of the approximation procedure is based on the following facts 
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Let F(co) be a. positive, even function on (-iTj-rr) (a "power spectrum") and 
let its Fourier coefficients be denoted by 


/ cos(ka))F(m)dw 


[these form the autocorrelation sequence for a stationary second order random 
sequence with power spectrum F(u)]. Let R(n) be the (n+l)x(n+l) Toeplitz 


matrix 



Then R(n) is positive definite for each n . Define [a^,p^(n) ,P|^(n)] 
to be the solution 




Then the following are true: 
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(i) 


det R{n) 
“n det ROi-l) 


> 0 


(ii) 

(iii) 

(iv) 


The 


1_ 

2if 


1 

2ir 


zeros of Pp(z) in the open unit circle, i.e. 


P„(x).= 0 => |X| < 1 

IT I IT 

/ cos{o)k)F (w) = ^ / cos(wk)F(w)dw = rj^ 

-tr -IT 


/V(e^‘^)P.(e"^‘^)F(aj)dw = 

-TT ^ 


0 , if 1f3 

I a. , if i=j 


for k=0,l,*»»,n . 

( Orthogonal polynomi al s ) 


(v) An equivalent version of the above is the following. Let T(n) be the 
(n+l)x(n+l) upper triangular matrix whose diagonal elements are ones, 
defined by 


T(n) 


ij 


0 , for 


i>j 

otherwise 


Then 

T^(n)R(n)T{n) = diag(a 

0 n 

This is a "triangular decomposition" of R(n) . 


(vi) is the minimum value of R(n)i|; , taken over all vectors 4' of the 


- form 


= 


X 

X 

X 


■ X 
1 
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(vii) The- sequence a ,a is nonincreasing and has the limit 

0 1 . • 

1 ' 

lim a = ^ F(o3)dw] 

- -TT 

, / 
(This is the geometric mean of the power spectrum F ). 


(viii) The solutions to equation (-1) may be computed recursively in n. . 
This is 

\ 

Levinson's algorithm : 


Initialization: 
Recursion : 


a ,= r 
0 ,0 


Pk("*lj ' Pk<"> - <^>Pn-n-k(") 

n 


In the above p^(n)=l for all n . 

Remarks : Levinson's algorithm dates from the paper: N. Levinson, "The Wiener 

rms error criterion in fitter design and prediction", J. Math. Phys. , 
Vol. 25, pp. 261-178, Jan. 1947. 


It is reprinted in the book Extrapolation, Interpolation, and 
Smoothing of Stationary Time Series by N'. Wiener, as an appendix. 


The result (vii) is due to Szego, and may be found in the book 
Toeplitz Forms and their Applications by Grenander and Szego. 
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For more references, proofs, and discussions see "The Use of Second- 
Order Intomation In Recursive Filter Approximation", by T.Mullis and R. 
Roberts, (Dept, of Electronics and Engineering, University of Colorado, 
Boulder). 



ThiS' is the discrete-time equivalent of the approximation problem. In 
other words, given a magnitude squared function F(w) , -ir<a)<7r , find a 
rational function 

H(z) - g ^ 

where 

' q( 2 ) = + + ..• +q^ 

p(z) = z" + + P„ 

m ^ n , the zeros of both p and q are in the open unit circle, 
such that 

The all -pole specialization is m=0 , so that q(z)=l . In this case a useful 
approximation is provided by Levinson's problem. If one sets 


g= 

pCz) = P^(z) 


I 


REPRODUCIBILnY OP OHE 
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Update Denominator 


Given F(o)} and q(z) , find an all -pole model 


a 


F((o) 


|p(ej“)|“ |q(ej“)|^ ’ 


set g = 


Update Numerator 

Given F(io) and p(z) , find an all -pole model 

^ - ■ . ■ set q = ^ 

■|q(eJ“)!^ F{a))ip(-eJ“)|^ ’ fST 

These two iteration modes are repeated until some degree of accuracy is • 
obtained. Note that one can increase the degrees m , n of the polynomials 

t 

between iterations. It is also possible to consider approximations of the form 

9 ^(zF-'-pgU) 

then the resulting approximation jH(e'^“)i^ is useful in that Fourier 
coefficients are matched, i.e. 

7T • 7T 

^ |H(e^‘^) |^cos(kw)d(o = ^ / F{u)cos{ko))da) 

-ir -7T 

for k=0,l,**»,n . Furthermore, the poles of H(z) are indeed in the open 
unit disc. 


Since low order all -pole approximations are not very efficient, it is 
advisable to consider the case m>0 , but how .can this be done? In RFA an 
iterative procedure is used. This procedure has a number of modifications, 
but basically .it goes as follows: 
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Initialize by setting q(z)'=p(z)=l , and update only one (factor) 
polynomial at a time, if it is necessary to have a suitably factored 
approximation.. 

Each iteration requires the computation -of a power spectrum, and a 
number of numerical integrations to obtain autocorrelation coefficients 
before Levinson's algorithm can be applied. It is here that most of the 
computation time is spent. For the computation of the power spectrums, it is 
necessary to compute |p(e’^‘^)j^ or |q(e^“)|^ . For this purpose the routine 
SPE ("Special Polynomial Evaluation") was written for the numerical integration 
the routines APA and JP were written. 

Description of SPE (P,N,CN,SN,X,Y) 

Here the inputs are 

/ 

P = (p^,p^.-**,pj 
CN = cosw 
SN = sinw 

and the outputs are x , y , where 

I 

X + iy = p(e^“) , where p(z) = + ••• 

^ The algorithm is fast in that it requires only about n. multiplications and ,2n 
additions. It is based on the following; 
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"0 ll rol To sine ■ 

A = , b = , C =. 

_ -1 2cos0 _ 1 _ _ -1 COS0 


then for k>0 


CA’^'^b = 


sin k0 
cos ke 


For the system 


i|^(0) = 


i|;(k+l) ■= Ai|^(k) + bp. 


„ 0 n ^ 

i^(n) = a" + Z a"^" bp 

1 k=i 


n+i-k 


Cip(n) = 


sin ne 

+ P, 

cos n0 ^ 


Im p(e^'®) 
Re p(ej®). 


sinC{n-l)05 

cos[(n-l)0] 


+ ♦ • • + p 


4) Rational Approximation on jai Axis 


The technique is to transform the problem to approximation on the unit 
circle via the bilinear transform 


= l±s s = (^1 
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Define the map o)^9 by s=jto , z=e'^® so that 


_ -sine _ 
1+COS0 


e = 2tan“^(^) 


(3a) 

(3b) 


The numerical integrations on (0,ir) in APA assume a uniform step 

size, so that 0 will take on the values , This generates 

0 0 

corresponding values of cd . Relation (3a) is thus used in TABCS to define 
the set of frequencies. The center frequency (i.e. for 0=^) will be 
(0=1 . RFA is set up ,so that is a power of 2, up to 2’=128 . The indexing 
is such that 

f 

(o(l) = 0 

(o(N /2+1) = 1 
0 

(o(N +1) = w 
0 

Note ; We can multiply all frequencies by the constant factor co defined 
' 0 

above, such that to(N /2+l)=to , which explain the name given 'to the 

0 0 

central ^frequency. 


The essential dynamic range', in deca 


, IS 


D = log 


10 


“(o(N ) ' 


1+cos 


IT 


= 21og 


10 - 


S'ln 


TT 


0 
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N 

16 

32 

64 

128 

256 

0 






D 

2.02 

2.62 

3.22 

3.82 

4.42 


The log frequencies are distributed symnetrically about the center frequency, 
with the density greatest at the center. In other words 

log[(o(k +j)] = -logCo){k -j)] 

0 0 

N 

where k +1 . 

0 2 

5) Transformation of Rational Functions According to Bilinear Map - Emphasis 
of Deemphasis 

Let F (w) _be a magnitude squared function on (o,“) . This general ies 
0 

the magnitude squared function 

F(0) = F (tan I) = F (w) 

0 0 

on . If F is the magnitude squared of a rational function with an 

0 

excess of e poles over zeros, then F(0) will have a zero of order 2e at 

IT . It is not wise to approximate F(0) in this form, since its inverse will 

have a pole of order 2e at it , and the numerator update will be ill -behaved. 

* 

Hence F(0) is multiplied by an, "emphasis", so that the result will have no 
poles or zeros on the interval {0 ,tt) . A natural choice of emphasis function 
is suggested by the bilinear transformation. Given a rational approximation 
for (0,ir) , we must be able to compute the corresponding rational approximation 




16 


for' (0,«>) . This is the reason for the following operator B , whose inverse 
is implemented in the program by the function subprograms BT and NORMLZ. A 
consideration of this operation will lead to a convenient emphasis 
deemphasis. 


For any monic (leading coefficient units) polynomial 7r(t) of degree 
k , define the monic polynomial 

.(BirXz) = i(|j]-)(z+l)'‘ 


(B'VKs) = n(^)(l-s) 


Suppose that F ('m)='];H(jw) ,• where 


«=> = 9 ^ 


and 


q is monic of degree m , zeros in open' RHP 


,p is monic of degree n , zeros in open RHP 


n -• m = e 


Then, 


■ H(s) = g— |-i R- (z+l 

. (Bj(z) 

For z=e^.® , jl+z |^=2(l+cos6) Therefore 


^m 


(Bp)(ej®) 
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Hence, we may as well use 


f — )® 

4+COS0^ 

for emphasis (see TABLE) and approximate 

|) 

on (0 ,tt) , note that F{0) has no poles or zeros on this interval. 

Suppose this is done, i.e. we have found an approximation 

H(z) = g 

p(z) 

where q is monic of degree m , with zeros in open unit circle and p is 
monic of degree n , with zeros in open unit circle, for which 

F(w) 


|ii(e^®)r . F(e) = 


C2(l+cos0)r 


How do we recover the approximation to F (w) which is stable and minimum’ phase? 


From the preceding approximation, we have 

F M = iH(j»)r » (ii+e^®r)® g“ 




Then 


H(s) 


(1+2)1'-"' g iM 
p(z) 

■1+S 


z =-^ 

^ 1=S 


_ , 2 xn-m ~ ^^1-s^ . r on-m ~ (B‘^q)(s) 


We may therefore take 



g = I (2)"“"' i I 

q(s) = (B'^q)(s) 
p(s) = {B‘'^p)(s)- • 
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Numerically, what is required in order to transform the jw axis problem 
to the unit circle problem is,, 

I 

(i) Emphasis, i.e. replace F (w) with 

0 . 

F (w.)- 

F(e.) = ° ^ 

^ [2(l+cose.)]" ^ 

(ii) After finding the approximation on (0 ,tt) , recover, g , q , p as above. 
In order to do this, the function subprogram BT was written. 

6) Implementati on 

a) Table (npq) 


For an excess e=NPQ ‘ we store a table in TAB of the values of 


u+cosw^ 


b) NORMIZ (FD,P,NPQ) 

✓ 

This subroutine transforms , F((o) in an emphaz.ised function FD(e-) with 
''-6€(0',Tr.). . 


REPR0DUCIBE;1TY OF THE 
ORIGINAL PAGE IS POOR 
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c) TABCS 


This routine computes a TABLE of Cosine and Sine. With step-size - 

^ n 

(n given) it computes cos(») , sin{-) and the frequencies 


W(-) 


SN(-) 
1+CN(*) • 


E. RFABND (TM,PNAP,PDAP.NNQ,NNP,GAMA) 


•Given the frequency response TH (input), a rational approximation 
GAMA*PNAP(»)/PDAP(‘) is calculated, with PNAP of degree NNQ , and PDAP 
of degree NNP . 

Note : RFABND begin to approximate with degrees (1/3) , then (2/4) ••• 

until a maximum error (EPS) of 1 % in the approximation achieved 
or until a maximum degree of 12 is reached- for the denominator PDAP^ . 

F- EXTLWB (NEXCESS) 

This routine tries to extend the Lower Bound of the overall transfer 

function using up to NEXCESS of poles over zeros. 

/ 

This is done by selecting suitable poles, which added to the upper bound 
still guarantee a correct behaviour in the time domain. 

Systematic cut and try is used here and so, starting from a first guess 
ISTART (index frequency) we converge toward the final value of ISTART , at‘ 
which one or several poles can be placed. 
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W('*) being the frequency-vector (see TABCS) then POLE=W(ISTART) . 

1) MUL (P.NP.POLE,PA,NPA) 

This routine calculates. PA(s)=P(s)*(s+POLE) where P is of degree NP 
and PA of degree NPA=NP+1 . 

2) KEEPS (BND.POLE.N) 

• f 

This routine updates the F.D.B. 

FDUM(-) = BND(-) * ( f 

P0LE^+W(*)^ 

FDUM is then compared to BNMN (BouNd Min) and BNMX (BouNd Max) which 

* i ' 

contain the lower and upper limit of the F.D.B. 

G. EXTHGB (NEXCESS) 

This routine EXTends the HiGher Bounds by placing appropriate zeros on 
the lowest bound TMIN and works in a similar way .to EXTLWB . 

H. FREQ 

RFA requires a frequency scale, given by TABCS . However Part II uses 
a logarithmic scale for the frequencies. The routine FREQ does this 
.transformation on the FREQuencies and makes the interpolation for the upper. 
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and lower bounds of the F.D.B., BNMN and BNMX and for the ratio R=BNMX/BNMN . 

Example : Suppose that NBR , the number of points/octave derived is 2 . 

Suppose k=2 , which means that the set of points considered in Part 
•II, is (1,1+22,1+2 x2^,1+3x 2S***) = (1,5,9,13,**‘) . In other words , 

j/ 

k=2 means that an interval of 2 =4 is chosen for the index 
frequency scale. 


Old system 
(TABCS) 


{ frequencies 
index- frequency 



New system 
(logarithmic) 


frequencies 

< 

index-frequency 



The logarithmic step is, ingeneral 



I. Other Subroutines Used 
1) INIT 

This is just an inialization of BNMN and BNMX , the lower and upper 
limits of the F.D.B. obtained with the second order model. 

k 

' The vectors FP(*) and FQ(») are also initialized - to 1 in order to 
be able to call the routine RFA . 
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At frequency OMEGA , the polynom P of degree N has for real part 
X and imaginary part Y, 

' P{jOME6A) = X + o'Y . 


WRODUCiBiLrrY OF -The 
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II. PART 2 
Philosophy 

The plant is i^'ead in (PLANT) and its range or templates (TMPLT) 
calculated over a set of frequencies. As an option, the boundary of the 
template can be plotted (P020) and the bounds on L can be computed and 
plotted (BNDG) . 

A. PLANT Model 


0 O 

1 2 


Fig. 2.0 

Between the nodes 1-2 , there can exist (Fig. 2.0): 

a) a known block of transfer function t(s) . This will be called a LINK ; 

b) an unknown gain k . In this case it will be called a KINK . 

The rules to describe the plant are the following: 

(i) KLINKSilO ; 

(ii) UKINKS<4 ; 

(iii) No feedback from final node; 

(iiii) No feedback into first node; 

(j) No trivial KINK , i.e., Uk<l ; 
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(jj) All integrations should be included in path from 1 to 2 ;'(The 
number of integrations is called NTYPE ); 

(jjj) *The only branch out from 1 is to 2 . 

1) Mathematical Preliminaries - 

Let us replace all integrations by a unit delay "1“ . 

Let IC(I) denote the length of the shortest path from node I to the . 
output. Thus IC(*) is a vector of the shortest path to the output, from 
different nodes. It is seen, then, that the shortest path from input to 
output is the excess e of poles over zeros of the plant P(s) . Let NPO 
and NQO denote respectively the degrees of the denominator and numerator 
of P{s) , so that e=NP0-NQ0 . 

Let lA be a matrix and IA(I,J) represent the length of node I to 
node J . 



Fig. 2.1 
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The minitnum path (Fig. 2.1) in 1 step only is represented by: 


lA 

1 


CO CO CO 

0 2 -1 

CO Co 00 CO 

00 0 OO 00 


The minimum path from node I to the output in 1 step is represented by the 
vector: 


IC = C = 

0 

Define addition (+) by 

lA (i,j) + lA (ij) = MinCIA (iJ),IA (ij)] 

12 12 

and multiplication (*) by 

lA (i,j) * lA (i,j) = lA (i,j) + lA (ij) 

12 12 

Then it can be shown that the set of matrices lA associated with the two 
laws of additive (+) and multiplicative (*) is a field. 



Thus, the minimum path in 2 steps is: 
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It can be proven that the minimum path in 2 steps to the output is: 

C = C + lA * C 
10 10 

and in general, the minimum path in N steps or less to the output is: 


Ct = * S-i 


vl 

It remains now to note that e=C^{l) when N is very large. 


2) Implementation 


“ in the program is replaced by the number 777 . The preceding formula 
(same notation) is used in the program to calculate e . 

The links are read in the matrix T , which is organized as follows: 


EEPEODUCIBILnY OP THE 
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link # 1 

’T(l), 

1 T(LINKS+1) i 

1 T(126*LINKS+1)“ 

link § 2 

• 

T(2) 

• 

! T(LINKS+2) 1 
1 . 1 

1 T(126*LINKS-tt2) 
1 ■ 

i , 

• 

link i LINKS 

• 

_T(LINKS) 

1 1 

1 * 1 
1 T(2*LINKS) , 

1 

i 

1 T(127*LINKS) 


frequencies W(2) W(3) W(128) 

T(i) is a complex number which represents the transfer function from node ~ 
k to A at the frequency W(n) . 

/ 

Then the links are read in and AK and BK are vectors of extremal - 
values for each KINK . 

B. RECQPG (NP.NQ,Q.P,G,X) 

For the plant condition X , this subroutine RECovers the polynoms Q 
and P of degree NQ and NP and the gain G , such that G*Q(S)/P(S) match 
exactly the. data contained in the previous matrix set. This problem is known 
as the Cauchy-interpolation problem and is presented under " Mathemati cal 
preliminaries “ below. 

1) Mathematical Preliminaries 


Suppose the data set {(A^,x^) ,(Xj^j,Xj^)} is given, where the {X^.} 
are distinct. 


tt(s) 


N 

tt(s-X^.) 



+ TT S 
1 


N-1 


+ 


+ TT 


N 


Define 


1 



Def i ne 






0 , k=0 

N xj"' 


k>0 


Then 

-k 

z h^s - = 
k=o 

7t(sJ 

and 


-1 

k= 

Construction 

U 

1—1 

1 1 

It 

O 

>] 


A = 

“ 0 

i 

I 



_-^N 

-IT 

1 _ 


Then 


b = 


A = 


SXkEk > 


1 


^^k'^ = FTxp- 


and 


cA*^”^b = ‘0 ' for k < N . 


Assume that the data satisfies 


where 


y^-p(X^) = q(X^.) for i=l,' 
p(s)' - s*^ + a '■s*^ ^ + ••• + a 


q(s) = q s"’ •+ + q ■ 


Note. : Distinguish between n and N . 
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0 , k=0 


a::' q{x.) 


, ■ k>0 


1 = 1 ' T 


E a s-'' - 


{9|,> = {0,T»*,0,q ,q +q 2X. 

^ 0 ^1 0 I 


k=N-m 


N xfVi 

"k ‘ iV 


Since y^p(X^)=q(X^. ) for each 1 , assuming that p(X.)?*0 for each i 


N x‘5-' q(X.) 

\ = .f ^ PTx^ 


- six) 1 _ ^ V x^"' 


ff-v\ _ q(x; 1 _ q(xj ^ 

TOFx ptiy^ 


^(x.) _i^ 

¥TI^ = 


Eh|^s'k = cq(A)[p{A)(sI-A)]‘*b 


n . N x|"*q(X,) 

j ir'(X^.) ° % 


This, implies the following 
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'h 

1 

*^n+i 


“®n" 


"S 1 

h 

2 

*^n+2 



= 

'S’a 

h 

3 

^n+3 


a 

1 


^3 




_1 








(3) 


Since gj^=0 for k<N-m , if N-m^n , i.e, N^n+m , then p(s) is determined 

from , {h . These numbers depend only on the data, and not on the 

{gj^} • Once p(s) is determined, the above equation can generate* the 

1 

sequence Which' can be used to generate q.(s)' in view of the identity 

. This is done as follows: 


qjs) 

^ sj 


l-'(N-m) X o .-(N-m+i) ^ X . o'N 
q^S ' , ^ + q^s ' ^ + ••• + q^s 


1 + ir s."^ + 


^ V 


-N 


therefore 


{0,0y« ,0,q ,q ,***,q ,0,0,0,.*»*} 

, ■ ^0 1-- m 

iindex=N-m 


= W.0,--,0,gj^_^*,g^_j^^,*..} * {1,7T 0,0,0, •••} 


i.e. 


^0 ■ ®N-m 


~ '^N-m+i ^i^N-m 

1 

\ “ ‘>N-m+2'^,‘'N-m+i * \8N-in 


etc. 


WRODUCIBiLrrY OP THE 
ORIGINAL PAGE IS POOR 



31 


Case 1 Ns2n‘ 


define 


then 


{y^-> = {(a^+jbj,*** ,(aj^,+jb^.),(a^-jb^),**‘,(a^,-jb^,)} 

n' (a.+jb.)(ju.)'^"^ 
h. = 2 2Re[ ~ . - n; ' 1 for k>0 

n' 

7r(s) = 7T (s^-hi)p 
1 

= 2 ji 0 j^Il(a)?-a)p CH is a deleted p ] 
n' 

h|^ = Z Re[Q^{bj^-ja^)(jw^) ] 


o 

It 


T~< 

11 

va 

k = 2 


k = 3 



Case 2 Same as case one except N=2n'+1 , add 0 to {X^} 

n' 

tt(s) = sn (s^+w?) = stt(s) 

1 

n' 2 

tt' ( 0) = (n OJ-) 

1 

= -2u2n(w?-w^) = jcOj^TT'-(ja)jj^) 


y 


=a 
0 0 
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n' a»+jb 




(n 10^) 


H=l 




k = 0 




01 


a 


^ = 1 “>A 

' 2- a^QjOi^ 

'' ' 3 -bjQjOi| 


.^oQc 


u 

Oloi.)^ 

- 2 

a 


, k>l 




k=i 


2) Implementation 


For the set of select frequencies W(*) , we compute P(j£o) . This is 
done through subroutine PJ . 


a) PJ (KW.X,PJW) 


Input : KW is index frequency 

X is the plant condition 

i 

Output: PJW <= PEjW(KW)] = P(jo)) 

In the frequency domain (Fig. 2.2), we can write: 

X = Gx + cu 


y = bx 
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Fig. 2.2 

Thus: s y - c(I-G)”^bu 

Here c==[l,'0] so the transfer function we are looking for is: 

PJW = [(I-G)"'b]^ 

where the index i refers to the first coefficient of the vector under 
parenthesis. Instead of inverting (I-G) , we employ "Grout reduction". 

Essentially this leads to the building of the matrix A<=CI-Gjb3 which 
is n+lxn where n+1 is the number of nodes. 

"Grout reduction" uses matrix A to compute 

PJW = C(I-G)"^b]^ as Cb+Gb+G^b+**0^ 

b)' J AKRON (KFLAG,NN) 

This routine builds up the Hankel matrix 
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^n+i *^2M+x 


used in equation (3). 
c) HANKEL (NUH,A) 


This subroutine solves recursively equation (3) and leads to the 
computation of the coefficients of the polynom P(s) of degree n . 

The sequence {gj^} then generates the numerator g(s) whose degree 
n-e=m is known, and the gain G . 

Note : The main reason for using Cauchy interpolation problem is to guarantee 
the minimum form of the plant P(s) . 

C. PATTERN (KW, KINKS) 

: index-frequency 

KINKS: number of independent parameters 

This routine computes the contour of the template of P at frequency- 
index KW and stores the necessary points. 
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1) Mathematical Preliminaries 


Using the bilinear theorem, we know that if one parameter k varies over 

(k ,k ) the plant transfer function describes in the complex plane, an arc of 
12 . . 

a circle. To represent an arc of circle in the complex plane- we use the mapping 

f 

R (E 


f(x) ^ with X e R 

a,b,c e £ 

If X e (0,1) we have an arc of circle. The relations between a , b , c • and 
f(0) , f(%) , f(l) are given by: 


f(0) = b 


and 


f(%) = 


f(l) = 


til 

f.i 

a+b 

c+1 


The relations between a , b , c and f(0) , f(l) , f'(0) are given by; 


f(0) = b 


f(l) = 


a-^b 

c+1 


and 


f'(0) = a - be 


[f (x)=-A±c_J 

(cx+1)^ 


Thus we can show that:. 


f(x) = f(0) + 
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with 

Notation : 



2) TOUR, (IN,M,N0) 


M is the number of KINKS 
M 

N0=2 is the number of vertices in the space of parameters. 
IN is a vector which summarizes the path taken. 



Fig. 2.3a 


Fig. 2.3b 
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All the corners (Fig. 2.3a) are reached by taking the path 1 (Fig. 2.3b) 
which can be described by the vector IN={1,2,1 ,3,1,2,1,3} (called a "tour") 
which indicates only which of the 3 independent parameters is changed. Thus 
IN describes a* "TOUR" . 

Two other tours are possible by single circular permmtation. 

In PATTERN the vector IDGE(.I,J) reveals whether the edge from I to J 
has already been examined CIDGE(I,J)>Oil , or not [IDGE(I ,J)=0] . In each 
edge, trois points are examined and the parametric linearization is found in 
ARC . 

3) ARC (.C,FP0,F0,FH,F1) 


Given FO , FH , FI . This subroutine computes FPO and C , in order to 
describe the arc-of circle f(z)=f(0)+f ' (0)^|:j:y with z e (0,1) . 

4) IVARC (C,FP0,F0,D,T) 

Given C , FPO , FO and D=f(l)-f(0) . This subroutine computes the 
curvature of the circle and depending on it, a multiple of 2 (points) are 
added to describe the arc of circle. 

IVARC can be 1, 2, 4, 8 , 

The eventual added points are set in the vector T-«-(t ,«**,t ) . 

I 8 
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+ 


Fig. 2.4 


5) KHULL (Z,N,X.I1.I2) 



Fig. 2.5 


Given some ordered vectors of a convex polygon S , Z(1)»»*Z(N) all 
complex, check if X complex belong to the convex HULL S . 


KHULL = i 


0 if X e S 

1 if 
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If WIULL=1 , II , 12 are the indices corresponding to the end values of the 
points to be next. (Example: here 11=1 , 12=3 ). 

6) IRT (X,Y,Z) 



Fig. 2.6 


Given X , Y , Z complex^ the value of IRT is returned according to the 

/ 

drawing. IRT is used by KHULL to detect if a point is inside or outside a 
convex hull* 

7) P020 (Z) 


This is a plotting routine for the templates in the Nyquist plane. 
P021(N,NPT,Z) plots on the Nichols' chart. 

8 ] POl (NQ.NP.Q.P.GAIN) 


This is a printing routine for the rational function GAIN*Q(S)/P(S) 


REPRODUCIBILrrY' OF THE 
OMGINAL PAGE IS POOR 
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9) DOT (-X.Y) 


t • t = DOT(x,y) 


D. TMPLT (KINKS) 


This routine first eliminates all the arcs which are interior to the 
convex hull returned' from the routine PATTERN . 


Then the contour of the template is stored in parametrized form. 



Fig. 2.7 

FO(KW) is the starting point of the template at frequency . NPT is the 
number of arc of circles which describe the template at KW . LOC is an adress 
such that FPO(LOC) contain the value of f'(0) for the arc no.l ; FP0(L0C+1) 
contains f'(0) for the arc no. 2 ; FP0(L0C+NPT) contains f'(0) for the 

arc no. NPT . 
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NTMP(KW) 


LOC and NPT are stored in arc register NTMP(KW) . They are obtained using 
the rule: 


LOC = NTMP(KW)/128 
NPT = NTMP(KW) - 128*L0C 

In the parametrization we do not use c , but store h=l/(l+c) in H(*) • 


Thus the first arc is described by: 

. FO(KW) , FPO(LOC) , H(LOC) 

It's final point is f^(l) and will be the initial point of the second arc 
described by: 


f^(l) , FP0(L0C+1) , H(L0C+1) 

and so on-. 

Note : Instead of storing the template of P , we store the template of 1/P . 
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III. PART 3 


Philosophy. 

We take a first guess for L (jo))' (FRSTG) and do 4 iterations of "local 
\ 2 

changes",. The "local" changes are done at each frequency and the correction 
then affects, all frequencies. This allow rapid initial ‘ball -park' convergence. 

"Global" changes are then made to improve the solution. "Global" refers 
to the fact that the corrections are made simultaneously at every frequency. 

1) Mathematical Prel iminari es 



Fig. 3.1 


Suppose that B(KW) is the bound on L(ju)) at the frequency. KW , and 

let M be the value being used at this point for L(KW) . The correction 

from point M- to point N is made by means of a lead-lag network, of the 

form t (s) = . 

cor' ' a s+b 
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Suppose we want the effect of the network to be maximum at the frequency 
W(KW) . Then we can let b=x . W(KW) and a=W(KW)/x and 


t^Q^tjW(KW)3l^ = x^ = y 


because 


t = x[W(KH)+xs3 

^cor^^^ [xW{KW)+s3 


We know that for minimum phase system , the phase and magnitude are 
related through the Bode-integrals. However, to use the Bode-integral we must 
evaluate singular integrals, difficult to do with precision. 


We use instead a technique first proposed by Bode {Bode, p. 343) whereby 

the given characteristic, is approximated by "semi-infinite segments". The 

phase corresponding to such a characteristic (break-point o) ) is given by: 

0 

+ 1f 0<x^..414 



log 



2 

TT 1+X^ 


if • .414 < x^ < 1 


where x = — . 
c w 

0 

Our procedure for global changes is as follows. Given L^^^(jcu) obtain 
after i iterations, we compute the Z corrections corresponding to the Z 
different frequencies of interest. Thus we will end with a set of Z 
magnitudes, from which the corresponding z phases will be obtained. The new 
set of z magnitudes and Z phases give new loop transmission , 


etc. 
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A. PLNMNL 


This routine calculates the Plant NoMiNoL for x=l , which correspond to 
the minimum value of all KINKS . The nominal jDlant is then stored, in 
PLNL{ ) for the magnitude squares, and in PLNFI( ) for the phase in radians. 


B.' FRSTG (GFJST2,1ST1,NNN) 

The first guess of L „„„ is takes as • Thus (L „ =GP„^„) , 

^ 2 nom s(s+a} 2 nom nom' 

r = k . 1 . 

^first guess s(s+a) P^q^(s) * 

GNVAR denotes variation in the high-frequency gain factor of the plant. 
ISTl is an index- frequency which is chosen such that 

' T 

R(ISTl) = > GNVAR but [R{ISTl-i)<GNVAR] 

'min I=IST1 

■Then a = W(ISTl) 


IST2 is an index-frequency calculated in Part I, and which essentially gives 
the first index-frequency for which L(jw) satisfies its bound. 

k is then chosen as the minimum gain -for which L satisfies its 

^ 2 nom 

bound on the interval (IST2,IST1) . 

The procedure now consists of finding which satisfies its bound 

on (IST2.,IST1) then to increase ISTl until than 

GNVAR (the high gain variation). The later will be the solution of our problem 
if we were satisfied with an excess of 2 poles over zeros. 

)F THE 
POOR 


BIPRODUCIBiLrrY « 
AMGINAL PAG® IS 
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Fig. 3,2 


C. BNDL (I$T2,I$T1,NNN) 

This routine computes and plots the first nine bounds on . 

D. XLDLAGl '(A,B,I»KL,NFCAG1,IST2>IST1) 

At frequency- index I e (IST2,IST1) , this routine computes the LeaD-LAG 
B 

Xs+B brings to its boundary. This is done by cut and try. 

RNGl (see below) computes and accordingly x is increased or 

decreased. 
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A and B are then returned to the current program. . 


Note: 


R = 


j2 

max 

T^. 

min 


1+L!^ 

'max 


ll+GP ^ 

' max 


+ -P 
P'max 


1+L 


min 


|1+GP| 


min 


G + i 


2 

min 


Thus in XLDLAGl we need to call -RNGl after translating the origin 
FO by G . 


PMAX <= [G + 


2 

max 


PMIN <= 


|G 


+ 


i|2 

P'min 


1) RNGl (FO,FPO,H,N,ZMIN,ZMAX, PMIN. PMAX) 


Given a boundary composed of N distinct arc of circles parametrized 
by F0,FP0(1) ,*•» ,FP0(N) ,H(1) • ,H(N) . This subroutine determines the 
points and , with the greatest and least distances from zero. 

2) CHKl (F,P,ZMIN,ZMAX, PMIN, PMAX) 

' P=|F|^ , F complex. 

I 

This subroutine does the following: 

If P ■< PMIli then P'=> PMIN 

If P> PMAX then . ' P => ZMIN^ 

If ' PMIN < P < PMAX , nothing happens. 
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3) MIDPNT (F1.D,H,FP0,F) 


If an arc of a circle is parametrized by FI , D , H , FPO this 
subroutine returns the number of special points contained in the arc. 


0 , 

MIDPNT = 1 , 

L2 , 


no special, point on the arc 
one special point is contained 
two special points are contained 



A and B ’ are special 
points 


4) KUAD (A,B,C.X1,X2,Y) 

Given the polynom p(s)=As^+Bs+C , this routine computes the roots of 
p(s)=0 and stores them in y . KUAD is then the number of roots Which lie 
in the interval (xl,x2) . 

5) ERROR (I,NFLAG,IST2.IST1) 

This routine computes the normalized Euclidienne distance of the point 
"to Its boundary at frequency I . Moreover, it stores the maximum 
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of all distances computed on (IST2,IST1) , with one db considered 
equivalent in 'distance' to 5 degrees. 

NFLAG=7 , indicates that the disturbance bound has been violated. 

Note : In XLDLAGl (except when called for local changes), the value of x 

is limited to (.95,1.05) (5% of changes allowed at each step) in 

order to avoid oscillations in the iteration. 

An additional refinement was found to be useful in the final stage of 

the iteration for eliminating oscillations. In Fig. 3.4 suppose the n-th 

iteration gives the points A , B for L(jo) ) , L(jo) ) respectively. Say 

1 2 

the correction corresponding to , is made for the new set of magnitudes 

(for the n+1 trial). This smaller magnitude gives a larger rate of decrease 

of magnitude near , driving L(jw^) in the direction of A' . No correction 

is persumably needed at because B is on the boundary. However, since the 

•magnitude of A' is now being used at w , if B is used at w then the 

1 2 

average rate of magnitude decrease near w is now less than before, i.e. 

2 

ma,g.(B-A')<mag.(B-A) , leading tO' phase lead at w , in the direction of B 

2 

To prevent this we should lower also the magnitude at o) , in the direction 

2 

of B' . This refinement is incorporated at present only along the high- 
frequency disturbance boundary. For example at C,D,**«,G in Fig. 3.4 the 
updates would normally be, say, .95, 1.0, .98, 1.02, 1.05 . Instead we use 
.95, .95, .93, '.95, 1.0 . 


REPRODUCIBrLm’ OF tHE 
' ORIGINAL PAGE, IS POOR 
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F. FIXG (A,B,L.ISTJST1,IST2) 


This subroutine is used for "global" changes to update the set of 
magnitudes squared. At index- frequency L e (IST2,IST1) the vector F is 
filled with . For the frequencies smaller than IST2., 

assumed to have a slope of 6dB/octave and for the frequencies greater than 
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ISTl , the last slope is extended indefinitely. 

6. SLOP {I1J2,NNN) 

Given a set of magnitude square F(*) over (11,12) , SL0PE(*) contains 

\ 

the relative change of slope over the previous interval at each frequency. 

H. PHASE (W,I1,I2,NNN) 

Using the previous results stored in SL0PE(*) , this subroutine computes 
the phase at frequency W , using the technique previously described. 

I. XGAIN (IST2,IST1,NM)' 


.This subroutine replaces G(*)<=XK*G(*) where XK is a pure gain, obtained 
such that as many points lie above the boundary below. 

J. Other Subroutines 


1) DSP (NPT,NY,X,Y,SYM,TITLE,KFLAG,Y1,Y2) 

This routine is a plotting routine. 
NPT : . number of points 

I 

Y1,Y2: Min and Max of the ordinates y 
NY : number of curves 
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X : vector of abcissal values 

Y : two dimensional vector of curves y[^=y{*,K) 

TITLE: vector of alphanumeric characters 

SYM : vector of symbols for up to 10 curves 

0 , if general 

KFLAG: 

1 , if we want to plot f(x,y)=0 

In the latter case of should be stored in DSPZ . 

2) DSP2 (A,Y1,Y2.X1,X2) 

Used when KFLAG=1 in DSP 

3) IFOR 

Used for scaling purpose. 

4) FIXBET (I) 

This subroutine calculates the maximum possible phase of L (jw) over 

2 

all frequencies and plant-conditions. 



IV. PART 4 


Philosophy 

At the end- of Part 3 we do not have as yet the required loop transmission 
L(jw) . We have for reasons of convenience a loop transmission (physically 
unrealiable) with has an excess of (2-^) poles over zeros, where 6j^ refers 
to the phase margin of L (see Fig. 4.1). 
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We therefore multiply this by H=- 


K(s-i-z)^ 


with g=l-orT > in order 

2 \ 2 . yu 


(s^+Zcoj^s+o)^)^ 

to obtain a realizable loop transmission with an excess of 5 poles over zeros. 

‘ 0 )** 

Obviously, we do not want to change the gain at low frequencies, K=~ . 

z^ 


We let w^=Az , and search over the three parameters z , A and t; , 
given 0j^ , in order to find an 'acceptable' characteristic. The search, over 
A and ? has been experimentally and charts of suitable A(0j^) and ?(0|v|) 
were prepared and stored in the program (see Fig. 4.2). Therefore, it is 
necessary only to search over z , given 9j^ , C(0[y|) and A(0n^) . This is 
done by cut and try. 


When z is found, L^{j“) is available in the form of numerical data, 
over a set of frequencies. A rational approximation is then made of 


L^(s) 


(s^+2z(0^S+£i)^)^ 


which is finite non-zero at s=0 and at s=<» . The reason for this modification 
is to retain- the double complex factor in the final rational representation of 
L(s) . 


A. CHKBND (IST,IST2,IST1) 

This routine uses the chart (Fig. 4.2), to compute ? and A . Then the 
value of z is adjusted by cut and try. 

The final results are stored in block ' [BETAlj . 
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Fig. 4.2 

Notations ; BETA <= 3 XLAMDA .<= -X 

A <= z QSI <= ? 

* 

B'., FREQl- .( 1ST) 

This subroutine changes our logarithmic scales for the frequencies into 

a trigonometric one as defined by TABCS . This is done in order to be able 

to take a rational approximation of L .On that new scale IGP is 

s 

\ 

interpolated and started in F(*) . 


C. INITl 


This subroutine initializes- FP and'.FQ , for use in RFA . 

% 


REPRODUCIBILrrY OF THE 
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D. RFAFNL 


This subroutine computes 


Ls(s) 


(s^+2zw^s+o)^) 


u: 


and approximates it by 6AMA*PNAP(S)/PDAP{S) , where PNAP is of degree NNQ 
and PDAP of degree NNP , Then is approximated and G(s) calculated. 


1) CNVLTW (A.NA.B.WB.C.KFLAG) 


This subroutine computes the product of two polynoms A and B of degree 
NA and NB and puts the result in C . 

E. PREFI (I$T,IST1) 

The data for the prefilter F is obtained from 




F = 


Min Max 


/ L \ / L \ 

’’^'l+L'Min'l+L^Max 


•^Mn^Max /,„ , 1, /TTT; 

G ^ P'Min ^ FMax 


Notations : 


FM <=■ F is stored in the block PREF 
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PMAX <= i G + I 


P'Max 


GMAG <= |G[ 


This subroutine, thus, prepares the data of |Fp . RFA is then used to 
obtain a rational approximation. 


PREFI prints out a summarizing table; 

OMEGA 

stand's for 

frequency 

DBl 

stands for 


DB 

stands for 


XLMIN 

stands for 

in. mag. square. - 

XLMAX 

stands for 

(L/1+L)n,g^ in mag. square 

DIST 

- A < 

•stands for 

tolerated disturbance in mag. squares 

PREF 

stands for 

value of the prefilter 

BNMN(KW)' 

is the 

actual min. bound of T at KW 

TMN 

is the 

tolerated 

TMX 

is the 

tolerated T,^g^^ 

BNMX(*) 

is the 

actual max. bound of T 


F. FREQ2 (NEXPF) 

This subroutine extends the data for the prefilter, assuming an excess of 
NEXPF poles over zeros for the prefilter. " 
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G. RFAPF (IST.ISTl) 


This routine approximates the data contained in FN~GAN!A*PNAP(S)/PDAP(S) , 
which is then stored back in the common block |FINPF[ , to be available for 
printing and plotting of the prefilter response. 



